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Abstract. We find a numerical self-consistent stellar model by finding the distribution function of a thin disk that satisfies 
simultaneously the Fokker-Planck and Poisson equations. The solution of the Fokker-Planck equation is found by a direct 
numerical solver using finite differences and a variation of Stone's method. The collision term in the Fokker-Planck equation is 
found using the local approximation and the Rosenbluth potentials. The resulting diffusion coefficients are explicitly evaluated 
using a Maxwellian distribution for the field stars. As a paradigmatic example, we apply the numerical formalism to find the 
distribution function of a Kuzmin-Toomre thin disk. This example is studied in some detail showing that the method applies to 
a large family of actual galaxies. 
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1. Introduction 

Few exact analytic stationary axisymmetric solutions for the 
collisionless Boltzmann equation are known in astrophysics. 
Most of the solutions correspond to spherical systems such 
as the polytropes and Plummer's models (Schuster 1883 
Plummer 119111 Chandrasekhar 19391, and lowered isother- 



mal models known as King models (Michie 1963; Michie & 
Bodenheimer 1963 King 1966 1 which describe well globular 
cluster data. Also, we can find self-consistent spherical models 



by using Eddington's formula (Eddington 1916). Furthermore, 
at least half of spiral galaxies and most of the elliptical galaxies 
appear to be triaxial stellar systems. Models for triaxial systems 
have been constructed numerically, usually end products of N- 
body collapse simulations (Aarseth & Binnev fT978l Miller & 
Smith ll979al b: van Albada fT9"S2l Wilkinson & James fT9"S2T l or 
using catalogs of numerically integrated orbits (Schwarzschild 
[T9791 Statler [T9871 Teuben fT987b : and analytically for spe- 
cial cases (Freeman l 19661 Vandervoort 1980). Moreover, other 
self-consistent triaxial models that satisfy the Stackel potentials 

Hunter 



were studied (de Zeeuw 1985, de Zeeuw et al. 1987 
& de Zeeuw 1992). Recently, a general solution of the Jeans 
equation for triaxial systems has been found (van de Ven et al. 
2003 ). But, in many astrophysical objects light is emitted by a 
thin stellar disk, so planar galaxies are of more than academic 
interest. The most important exact analytical solutions of the 
collisionless Boltzmann equation for thin disks are the ones due 
to Mestel dl963> and Kalnajs J1972i . Also, Miyamoto (119711 ) 
found an approximated analytical solution for the Toomre disk 
using a series method, and Ng ( 1967) found a self-consistent 
model for a thin stellar system, from a given distribution func- 



tion, performing numerical calculations. The importance of 
having stationary solutions, for stellar objects or other systems, 
of the Boltzmann or Fokker-Planck equations is that they are 
a starting point for doing dynamical evolutions, and to study, 
analytically or numerically, the stability of the system by per- 
forming different kinds of perturbations. These studies have 
been done over the years by many authors, as examples we 
can mention (Cohn 1979 Watanabe et al. 119811 Nishida et 
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al. lT9g51 Yepes & Domfnguez-Tenreiro lT9921 Takahashi [T995l 
Theuns [T996l Takahashi et al. [T997l Joshi et al. l200Tl Ashurov 
2004). In most of the thin disk models, we only have as in- 
formation their Newtonian gravitational potentials and surface 
mass density, knowing its distribution function we can obtain 
the thermodynamic variables allowing a full understanding of 
the physics of each individual model. Thus, the solution of the 
Fokker-Planck equation for each particular disk is needed to 
determine its main physical properties. Despite of the impor- 
tance of the Kalnajs and Mestel thin disks, we have to keep 
in mind that they do not model a realistic physical situation 
because they exhibit odd behaviors within the disk, i.e. in the 
Kalnajs disks, the disttibution function / goes to infinity for 
certain values of the energy density and angular momentum; 
and in the Mestel disks, the gravitational potential 4> diverge at 
the center of the disk. 

The Fokker-Planck equation is also known as the Fokker- 
Planck approximation because it truncates the BBGKY 
(Bogoliubov, Born, Green, Kirkwood, and Yvon) hierarchy of 
kinetic equations, at its lowest order, by assuming that correla- 
tion between particles only plays a role as a sequence of uncor- 
rected two-body encounters. It is worth noticing that the only 
"approximation" made in the Fokker-Planck equation comes 
from the model adopted for collisions and, in fact, the Fokker- 
Planck equation with a general collision term can be derived 
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from first principles and no ad hoc suppositions are needed. 
The main goal of our work is to find and study a self-consistent 
model of a stellar thin disk, i.e. a distribution that satisfies 
simultaneously the Fokker-Planck and the Poisson equation, 
the meaning of simultaneously will become clear later in the 
manuscript. In writing the Fokker-Planck equation, we consid- 
ered in the collision term the local approximation and used the 



diffusion coefficients found by Rosenbluth et al. ( 1957 1. 

The computational problems in solving the Fokker-Planck 
equation without using any statistical (Monte Carlo method) 
and average (moment equations) are briefly discussed in this 
paper. The Fokker-Planck equation is solved by a direct nu- 
merical solver using finite differences and a variation of Stone's 
method (Stone 19681 Ujevic & Letelier 2005 1 to include mixed 
derivatives. This allows high resolution solutions with lower 
computational cost (Ujevic & Letelier 2005 1. Due to nature of 
the problem, either analytical or numerical solutions are very 
difficult to obtain. In this respect, it is illustrative to cite the 
words of Binney and Tremaine ( 1987 page 245): Finding the 
particular function of three variables that describes any given 
galaxy is no simple matter. In fact, this task has proved so 
daunting that only in the last few years, three-quarters of a 
century after Jeans' s (Jeans IT915I paper posed the problem, 
has the serious quest for the distribution function of even our 
own Galaxy got underway. As far as we know, there are no ana- 
lytical or numerical solution of the Fokker-Planck equation for 
thin disks. 

The article is organized as follows. In Section II we present 
the Fokker-Planck equation for a thin disk with the collision 
term given by the local approximation and the diffusion coeffi- 
cients found by Rosenbluth et al. \\951\ . In Section III, we ex- 
plain briefly the computational codes used to solve the Fokker- 
Planck equation. In Section IV, we present the properties of the 
Kuzmin-Toomre thin disk used as an example for the numerical 
method. The numerical results for the velocity and spatial dis- 
tribution functions for galaxies based in this model are showed. 
Finally, in Section V, we summarized our results. 

2. Fokker-Planck equation for a thin disk 

The general problem for the evolution of a system made with 
stars of mass m using the Fokker-Planck equation consists in 
the determination of a distribution function / and the gravi- 
tational potential O that solve the Fokker-Planck and Poisson 
equations, say 



^- + v-v/-vo-v„/ = r[/], 



V 2 4> = AnGm 



\f d \ 



(1) 
(2) 



where v represents the velocity of the stars, V is the usual spa- 
tial gradient, V r is the velocity gradient (derivations are done 
with respect to the velocities), G is the gravitational constant 
and the symbol F[/] represents the Fokker-Planck collision 
term. The solution of the Fokker-Planck and Poisson equations 
is found from an iterative procedure in an short time interval 
with respect to the time scale in which O varies. The method 



of evolution has two main steps: (i) the distribution function is 
first advanced in time using equation Q with the potential held 
fixed, and (ii) the potential is then adjusted so that the Poisson 
equation is again satisfied. Some systems, as a final state of the 
evolution, may reach a stationary regime. Then, the temporal 
part in Q is equal zero and the distribution function is time 
independent. Hence, in this case, the distribution function must 
satisfy simultaneously the stationary Fokker-Planck equation 



vv/-vo-v l ./ = r[/] > 



(3) 



together with equation 0. Numerically we have that, in a non- 
stationary regime, a gravitational potential at time t creates 
through the Fokker-Planck equation a distribution function at a 
later time t+At. In our case, a stationary regime, we want to find 
a gravitational potential and a distribution function at the same 
time t. This is what we mean when we say that the solution 
(O, /) satisfies simultaneously the Fokker-Planck and Poisson 
equation. Note that we can start with a given gravitational po- 
tential at time t and make the evolution of the system by first 
finding its distribution function at a time t + At, and later, the 
new gravitational potential, via equation Q, at this new time 
t + At, and hope that the successive iterations of this procedure 
converge. In the case of convergence, we obtain also a station- 
ary solution but this solution might not have any physical sig- 
nificance, i.e. the final gravitational potential together with the 
distribution function do not model any astrophysical system. 
For that reason, it is important to find the simultaneous station- 
ary solution of the Fokker-Planck equation for a given gravita- 
tional potential that can describe a real system. The procedure 
to find this solution is explained in the next section. Our main 
goal is to obtain a physical consistent model by solving the set 
of equations J2I3> without using any approximations, like the 
ones found in the orbit averaged approximation, Monte Carlo 
method, or the Jeans equations (Binney & Tremaine 1987), i.e. 
a direct numerical solution. 

To find the stationary Fokker-Planck equation valid for thin 
disks, the non-stationary case is straightforward, we start from 
the general case in three dimensions and later perform the pro- 
jection to the thin disk plane. The collision term in the local 
approximation has the form 



r[f] = -T J —[f(x,v)D(Avd] 
2 f— { dvjdvj 

1,7=1 



(4) 



where the functions £>(Av,) and D(Av,Avj) are known as 
the diffusion coefficients. These coefficients were calculated 
(Rosenbluth et al. 11957> assuming and inverse square force 
for the interaction, and also, that each stellar encounter involve 
only a single pair of stars and are independent of all others. 
They are simplified if the field stars distribution function is 
a Maxwellian distribution (Binney & Tremaine ll987> . Recall 
that the diffusion coefficients are calculated considering a test 
star of mass m moving through an infinite homogeneous sea 
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of field stars of mass m a who has mean velocity equal to zero. 
We know from stabilities analysis that this gravitating system 
can not be in static equilibrium and to overcome this problem 
we invoke the Jeans swindle. In order to calculate the diffu- 
sion coefficients in the next sections, we are going to set the 
dispersion of the field stars equal to the typical velocity of the 
problem considered. 

The diffusion coefficients described above are valid for a 
general three dimensional case. To obtain the Fokker-Planck 
equation for thin disks, we perform a projection on the z — 
and v, = planes using Dirac deltas. Doing the replacement 

f{x, y, z, v x , v y , v z ) -» fix, y, z, v x , v y , v z )6iz)6(v z ), (5) 

in Q, and integrating in the z and v z variables, we arrive at 
a stationary Fokker-Planck equation in two dimensions, see 
Appendix |X| The indexes (i,j) of the collision term @ now 
take only the values 1 and 2. The important thing about this 
procedure is that the form of the Fokker-Planck equation and 
the diffusion coefficients are preserved but projected into the 
z — and v z = planes. 

3. Numerical method 

The solution of the Fokker-Planck equation in a stationary 
regime is not an easy task. The problems involve in the so- 
lution are both numerical and physical. Numerical because in 
the three dimensional case, the Fokker-Planck equation has six 
independent variables, three space coordinates (x) and three ve- 
locity coordinates (v). In the two dimensional case, the disk 
case, a simplification exist because the total number of vari- 
ables are four. In either cases, in a finite difference scheme, 
the large number of grid nodes needed for the computation of 
the solution becomes a data storage problem. In a three dimen- 
sional problem, if we divided each of the six variables intervals 
of the distribution function in nine parts (ten nodes), we will 
have a grid with 10 6 nodes. This corresponds to a same num- 
ber of linear equations to solve. In a simple numerical method 
we have to store and solve a matrix with 10 12 elements. For 
two dimensional disks, the main matrix will have 10 8 elements. 
So, our problem is slightly simpler than the three dimensional 
case. With 10 grid nodes per variable only very simple geome- 
tries can be described. Furthermore, the large number of matrix 
elements brings us another computational problem, the slow- 
ness of the codes. Most of the elements of the these matrices 
are zero and this motivates the search for alternative and faster 
methods to solve the linear system using only the non null data. 
We overcome these problems performing a discretization of the 
Fokker-Planck equation using a finite difference thirteen-point 
molecule method and then solving the resulting discretization 
matrix with a modification of the method proposed by Stone 
(1968 1. The original method of Stone was modified to han- 
dle two (four variables) and three (six variables) dimensions, 
and the extra diagonals that appear when mixed derivatives 
are used (Ujevic & Letelier 20Q3. Also, it can handle non- 
stationary problems, and curve boundary conditions with a lit- 
tle modification of the code. Remember that the collision term 
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in the Fokker-Planck equation has mixed derivatives. Physical 
because the search for a boundary condition is not always clear. 
In general, the analytic potentials that we have are only geomet- 
rical, i.e they have been developed without using any observa- 
tional data (only an image method solution for a point mass). 
From this fact, it is not obvious that these potentials could rep- 
resent some real galaxies. 

The surface mass density of the disk is obtained integrating 
the right hand side of (|2ji, say 

Xoo 
fd\. (6) 
oo 

Due to the fact that the integral limits of are from -oo to oo, 
we must define first a criteria to establish the velocity domains 
for the calculation of the distribution function. In general, we 
are going to use the value of the escape velocity v esc = y2\<S>\ 
at the center of the disk to establish the upper and lower limits 
of the integration. 

When Cartesian coordinates are used, the spatial disk part 
of the distribution function will be delimited by an irregular 
boundary (circumference). By using cylindrical coordinates we 
can transform the irregular domain into a more simple square 
domain, but in this case we need more boundary conditions and 
assumptions for the distribution function. For example, if we 
are using a Dirichlet boundary conditions, we need the value 
of the distribution function on the side of the square that cor- 
responds to the coordinate R = ^x 2 + y 2 = (the center of 
the thin disk), a value that is not always known. To avoid these 
difficulties, we prefer to use Cartesian coordinates in the spa- 
tial domain and treat the problem of the irregular boundary, the 
code presented in (Ujevic & Letelier 2005 ) allows this option. 
In compensation, we only need the value of / in the exterior 
radius of the disk, i.e. only one expression for the boundary 
condition. 

The integration program was tested by integrating expo- 
nential functions of the form exp[(— v\ - v 2 )/2cr 2 ] for differ- 
ent values of <x. The integral is performed using Romberg's in- 
tegration scheme with Richardson style extrapolation (Gerald 
& Wheatley 1994 1. The Fokker-Planck code was tested com- 
paring the results obtained with the modified Stone method 
with other methods, as for example, the explicit Euler method 
and Gauss elimination, obtaining full agreement. This was per- 
formed for a small number of nodes because for larger grids the 
usual methods are not convenient. All these tests assure us that 
our codes are working properly. The LU decomposition used 
in the modified Stone method is presented in the AppendixlBl 

The general procedure to find the distribution function of 
a thin disk is as follows: first we start with a given stationary 
disk potential which is inserted into the Fokker-Planck equa- 
tion. Then, we find numerically a distribution function in phase 
space associated to the gravitational potential and the bound- 
ary condition imposed. Later, we perform the integration (|6j 
using the recently calculated distribution function to find the 
surface mass density of the problem. This surface mass density 
has to be equal to the surface mass density of the given grav- 
itational potential if we want to have stationary self-consistent 
solutions between the Fokker-Planck and Poisson equations. If 
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Fig. 1. Nodes in the spatial domain to be used in the two di- 
mensional figure of the surface mass density. The full circles 
represents the nodes inside the disk used for the calculation. 
The empty circles are necessary for maintaining an order in the 
grid but do not enter into the solution. After finding the surface 
mass density, the spatial part of the solution of the distribution 
function can be divided into vectors (represented like arrows). 
This procedure help us to present the surface mass density re- 
sults in a simple way using a two dimensional figure. 

the surface mass density obtained after the integration of is 
not equal to the one given a priori then it is an indication that 
the boundary condition used for the distribution function is not 
compatible with the problem studied. Our goal is to find the 
self-consistent solution of the problem by adjusting the param- 
eters involve in the gravitational potential, galaxy model and 
boundary condition of the Fokker-Planck equation. Note, that 
in this procedure we do not solve the Laplace equation. This 
was done only for didactic reasons and because there exist in 
the literature a huge variety of methods that solve this equa- 
tion. Actually, to obtain a self-consistent model, we only need 
the surface density of the thin disk. With this information, we 
find the gravitational potential via Laplace equation, and later 
we apply the above procedure. 

Now, we need a boundary condition for /. An appropriate 
boundary condition can be found if we recall that the veloc- 
ity distribution in a star system is determined by two compet- 
ing processes: relaxation through stellar encounters that drives 
the distribution toward a Maxwellian form, but stars beyond 
the finite escape velocity continually disappear from the sys- 
tem. Since escaping stars traverse the system in a small frac- 
tion of a relaxation time, the resultant velocity distribution 
should be zero beyond the escape velocity but otherwise nearly 
Maxwellian. Also, most of the stars in thin galaxies are cen- 
trally concentrated. So, for a stationary system, these consider- 
ations lead us to use as an appropriate boundary condition for 
our disk the function 

fBoun <* exp[(-x 2 - y 2 )/2cr 2 s ] exp[-(v - v typ ) 2 /2cr 2 ], (7) 

where cr s and cr v are the spatial and velocity dispersion respec- 
tively, and v = ^jv 2 + v 2 . The velocity dispersion <x,, is a mea- 
sure of the deviation from the star typical velocity of the galaxy. 



The spatial dispersion <r s is a measure of how centrally con- 
centrated these stars will be. The velocity part of the boundary 
condition is centered at the typical velocity of the system. 

4. Kuzmin-Toomre Galactic Model 

The code used, and the procedure presented in this manuscript, 
can be applied, in principle, to find the distribution function for 
any thin disk potentials. But, as an example, we are going to 
used the Kuzmin-Toomre family of thin disks (Kuzmin 1956 



Toomre 1963 Binney & Tremaine 1987 1 which, in cylindrical 
coordinates, has the Newtonian gravitational potential, 



GM 



y/R 2 + (a + \z\) 2 



(8) 



When z < 0, Q>tk coincides with the potential of a point mass 
M located at the point (R,z) - (0, a); and when z > 0, fl>rjr is 
exactly the potential of a point mass at (0, -a), hereafter a is de- 
noted as the cut parameter. The cut parameter modifies physical 
properties of the disk, like energy density, pressures and sound 
velocity profiles (Vogt & Letelier 2003 1. By applying Gauss's 
theorem to a flat volume that contains a small portion of the 
plane z — 0, we conclude that the potential (jSJi is generated by 
the surface density 



aM 



2n(R 2 +a 2 )i' 



(9) 



The potential describes a disk which extends from R = to 
R - co. In this work we shall consider a disk with finite radius 
with the same type of potential. 

We start with a galaxy model of radius 40 kpc containing 
10 12 stars of solar masses M . The typical streaming veloc- 
ity considered is of order 200 km s _1 with velocity dispersion 
cr v = 40 km s _1 . Hereafter, we set the masses of the field and 
test stars equal (m = m a = M Q ), The mass density p that appears 
in the diffusion coefficients is calculated spreading the total 
mass of the system in a thin disk, which radius is given by the 
model. To solve the Fokker-Planck equation for the Kuzmin- 
Toomre disk type potential we use a total of 810000 nodes, 
30 4 , distributed along the computational domain. All the nodes 
are necessary to maintain a computational order that is needed 
by the code, but actually we do not used all of these nodes in 
the calculation because some of them lie outside the thin disk 
region, see Fig.^ Recalling that, we have to check the results 
of our numerical solution of the distribution function indirectly, 
we compare the surface mass density obtained from the integral 
with the exact solution l|9}. This will give us an estimate of 
how far we are from the exact distribution function. In Fig.^ 
we present the spatial part of a disk with a given number of 
nodes used for the computation. 

For every internal node of the spatial part (full circles in 
Fig.^ we perform the integration of the distribution function 
obtaining in that way the surface mass density of the disk. 
Note that the surface density is only known in these internal 
nodes. This surface mass density can be plotted in a conven- 
tional three dimensional plot (x,y,Y). But, this makes the vi- 
sual comparison with the exact result Y,tk cumbersome. For 
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Fig. 2. Top: Two dimensional plot of the numerical solution for the mass surface density of the Kuzmin-Toomre like disk with 
cut parameter a — 40 kpc, and <t s = 35 kpc. The solid line is the exact solution obtained from the mass surface density of the 
Kuzmin-Toomre disk type. The dotted line represents the mass surface density obtained after the integration of the numerical 
solution, the full circles are the values of the numerical calculations. The arrows from Fig.^ are plotted one along the others. 
Each point represents a position {x,y} on the surface of the disk. Bottom: the same of Top but with cut parameter a = 62.7 kpc 
and cr s = 40 kpc. Adjusting the values of the cut parameter a and the spatial dispersion cr s we can approach the numerical mass 
surface density profile to the exact solution. 



ease in plotting and comparison, we shall plot the surface mass 
density in a two dimensional plot ({x,y},£). The nodes {x,y} 
are chosen following the direction of the arrows, starting from 
the left start arrow. Once finished, we moved to the next right 



arrow and plot the {x, y] nodes of this arrow along the nodes of 
the previous arrow. We repeat this procedure until we reach the 
end arrow. The nodes that are outside the disk (empty circles) 
have surface mass density equal to zero and will be omitted for 
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Distribution Function (a) 




Distribution Function (b) 




Fig. 3. Numerical solution of the velocity distribution of the 
distribution function for a Kuzmin-Toomre disk type when a = 
62.7 kpc and cr s = 40 kpc. The velocity distribution functions 
correspond to radius r = 39.26 (a), 20.37 (b), 1.95 (c) kpc. 



the plot. We use this procedure for the numerical solution as 
well as for the exact mass surface density solution. In the later 
case, we take the exact values corresponding to the grid nodes 
used in the numerical solution. We show the results in a two di- 
mensional figure, instead of a three dimensional one, because 
we think that the value difference of the mass surface density 
between the grid points in the numerical and exact solutions 
can be noted easily. The spatial difference between arrows is 



Distribution Function (a) 




Distribution Function (b) 




Distribution Function (c) 




Fig. 4. Numerical solution of the spatial distribution of the 
distribution function for a Kuzmin-Toomre disk type when 
a = 62.7 kpc and cr s = 40 kpc. The spatial distribution func- 
tions correspond to velocities v = ^jv 2 x + = 341.6 (a), 207.3 
(b), 27.3 (c)kms- 1 . 



the spatial discretization Ax used in the solution of the Fokker- 
Planck equation. 

In the upper graph of Fig. [2 we present the numerical re- 
sult of the mass surface density (|6} with parameters cr s - 35 
kpc, cr v = 40 km s _1 and a — 40 kpc (dotted line), compared 
with the exact density function (solid line) obtain from the 
data of the Kuzmin-Toomre density function (|9jl- In this case 
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Vesc ~ 460 km s and the integration have been performed in 
the intervals v x e [-350, 350] km s" 1 and v y e [-350, 350] km 
s _1 . By varying the cut parameter a and the spatial dispersion 
cr s , we can lower or raise the height of the ridges in the com- 
puted curves, note that this modifications also leads to changes 
in the escape velocity and integration intervals. In the lower 
graph of Fig. |2] we show the case with parameters cr s = 40 
kpc, o\, = 40 km s and a = 62.7 kpc. For these parameters, 
v esc * 370, v x e [-250,250] km s" 1 and v y e [-250,250] km 
s -1 . The maximum relative error between the computed val- 
ues and the exact values is in this case less than 1%. This is a 
good result because we are trying to match a real finite physi- 
cal galaxy to an infinite geometric potential that was developed 
without using any observational data. In solving the Fokker- 
Planck equation we must take into account that the final nu- 
merical result of the distribution function has to be positive. 
Numerically, this means that all the values of the distribution 
function at the grid nodes have to be positive. Numerical sim- 
ulations show us that this leads to restrictions in the velocity 
dispersion and velocity intervals. For example, if we lowered 
the values for the velocity dispersion to values near zero, the 
distribution profile becomes very sharp and then some of the 
values of the distribution function are negatives. The value of 
the velocity dispersion depends on the galaxy model under con- 
sideration. In particular, for the solution found above, we can 
not say that the distribution function is unique because it de- 
pends on the relative error allowed for the numerical solution. 
For example, if we allowed a relative error of 1.5%, the param- 
eters a and <r s can take values between [62.4 kpc, 63 kpc] and 
[39.9 kpc, 40.1 kpc] respectively; and cr v can not take values 
lower than 25 km s (because some values of the distribution 
function are negatives) but might have higher values. In gen- 
eral, for the gravitational potential 0, we were able to model 
galaxies with positive distribution functions for different pa- 
rameters values, i.e. galaxies with total number of stars ranging 
from 10 10 to 10 12 , typical velocities from 100 km s to 300 
km s _1 , velocity dispersions from 20 km s _1 to 80 km s _I and 
galaxies radius from 10 kpc to 80 kpc. So, in principle, we can 
find stationary solutions for many typical galaxies. 

Note that the potential we are using is from an infinite 
thin disk, and that we have made a cut at an usual radius for 
galaxies. So, the results we obtain are remarkable because we 
modeled a stationary galaxy with physical parameters from 
a infinite-cut disk which do not necessarily describes a real 
galaxy. Better agreement can be found if we enlarge the value 
of the spatial dispersion above the fixed radius of the galaxy. 
This behavior is equal for every size of the galaxy we choose, a 
numerical indication that the galaxy is actually bigger than our 
chosen radius. This is not surprising because the gravitational 
potential used in the Fokker-Planck equation is originally from 
an infinite thin disk. 

All the calculations have been done without actually know- 
ing the proportionality constant appearing in 0. To determine 
the proportionality constant, we take one point of the initial 
data for the problem, the mass surface density 1,tk 0, and the 
numerical mass surface solution at the same position. Later, the 
numerical solution is scaled with this constant. 



The numerical solution of the distribution function obtained 
in the cases mentioned before can only be seen by a collec- 
tion of three dimensional figures. For example, if we want the 
velocity distribution for the spatial part, we will have one fig- 
ure for each pair (x,y). In Fig.|5J we show the velocity distri- 
bution function at three different radius for a Kuzmin-Toomre 
thin disk with parameters <x s = 40 kpc, cr,, = 40 km s and 
a = 62.7 kpc. From these figures we note that, in order to have 
a stationary galaxy, the velocity distribution function have to be 
centrally concentrated. In all the cases, the population of stars 
with velocities that are near the escape velocity of the system 
is almost equal to zero. Also, in Fig.|4]we present three spatial 
distribution functions for the model of Fig. [3] We note from 
these figures that in general the population of stars increases 
with lower velocities, and that for each particular velocity the 
concentration of stars are bigger near the center of the galaxy 
than in the edge of the galaxy. Note that the distribution func- 
tion can have other symmetries besides the axially symmetry 
presented in the gravitational potential ©and surface density 
due to the form of equation 0. 

Note that the distribution function obtained from has 
to be positive and also, after the integration 0, the surface 
density profile obtained has to match the initially given sur- 
face density 0. These facts impose restrictions on the bound- 
ary conditions and for that reason the boundary condition can 
not be set arbitrarily. The boundary condition used for the 
numerical calculation satisfies all the above requirements. We 
have tried other boundary conditions by simple modifying the 
exponents of equation and with other different functions 
without obtaining good results. 

5. Conclusions 

We describe a method to find numerically a distribution 
function that simultaneously satisfies the Fokker-Planck and 
Poisson equations providing a self-consistent galaxy model for 
thin disks. This solution have been found by solving the orig- 
inal Fokker-Planck equation avoiding simplification using a 
previously developed code based on a direct finite differences 
method. 

As an example, we find distributions function of the family 
of Kuzmin-Toomre thin disks. As far as we test, we are able to 
model many self-consistent stationary galaxies using the most 
common physical parameters known for galaxies. This is ac- 
complished by adjusting the parameters found in the gravita- 
tional potential, galaxy model and boundary condition. In prin- 
ciple, we can find consistent astrophysical model for other thin 
disk potentials. Once the solution is found, the stability analy- 
sis of the models can be achieved performing perturbations to 
the stationary solutions previously found and using the same 
code presented in Ujevic & Letelier (20051 and performing a 
Crank-Nicolson discretization in time. We also found that 
is an appropriate boundary condition for the case of a stationary 
Kuzmin-Toomre stationary disk. In forthcoming articles this 
condition will be checked for other thin disk models as well 
as three dimensional disk models. 

In this manuscript we presented the results for a case in 
which the gravitational potential and the surface density is 
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known a priori. This was done only for didactic purposes be- 
cause, in general, the method presented here can be applied 
to more general problems in which we only know the sur- 
face mass density of the disk. With this condition, the grav- 
itational potential can be obtained, numerically most of the 
time, through the Poisson equation 0, and then the distribu- 
tion function / through the Fokker-Planck equation. Later with 
/, the complete surface density is found via the integral (|6). 
Then, the problem is completely solved by matching the sur- 
face mass density obtained via the integral (|6} with the one 
found from the Poisson equation. 

For future applications, it will be very important to test this 
code and boundary condition in other thin disk models and 
three dimensional structures in order to calculate their distri- 
bution functions. Three dimensional models, can be obtained 
following the procedure of this manuscript together with the 
code in Ujevic & Letelier (2005). In the last reference, a three 
dimensional calculation of a test Fokker-Planck equation is pre- 
sented. The construction of a solution for Plummer and King 
models, thick disks, and the stability of thin disks and other 
three dimensional structures are currently under investigation 
by the authors. A stability analysis of several thin disks using 
a perturbation on its energy-momentum tensor in the context 
of general relativity can be found in Ujevic & Letelier (2004i, 
another stability analysis of related structures using Rayleigh 
criteria of stability can be found in Letelier (2003 1. 



Acknowledgements. M.U. and P.S.L. 
support and P.S.L. also thanks CNPq. 



thanks FAPESP for financial 



Appendix A: Projection into the z-v z plane 

In this section, we develop some of the calculations that we en- 
counter in the substitution of the modified distribution function 
(JSJl into the stationary Fokker-Planck equation 0, 



II dz 



[f6(z)S(v z )]dzdv z = J 
= 0, 



?f$(z) + f6(z) 
az 



dz 



II 



df 
dz 



df 
dz 



z=0 



[fD(Av z )5(z)6(v z )]dzdv z 



li 



JL[fD(Av z )]5(v z ) + fD(Av z )5(v z ) } dv : 
ov z 



= —[fD(Av z )] 

ov. 



v,=0 



- —[/£>( Av,)] 
ov 2 



= 0, 



v.=0 



where we used the well known Delta function properties: 

X+oo p+oo 
-co g(x)S(x-a)dx = g(a) and J ^ g(x)6(x-a)dx = -g(a). 

The remaining projections can be calculated in a similar way. 
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Appendix B: Two dimensional UJ decomposition 

The LU decomposition for the two dimensional version of the 
modified Stone method presented in Ujevic & Letelier (2005 1 
is, 



where X p is the element p of the diagonal Y in matrix X. These 
relations have to be calculated in the order specified above. 
The notation used for the thirteen point molecule is describe 
in Tablelim 
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Table B.l. Relations and nomenclature between the matrix form and the one dimensional storage index at node p of the thirteen 
terms used in the discretization of the Fokker-Planck equation. 



Matrix Form 


Abbreviation 


Name 


Position From Node p 




Basic Diagonals 




f{Uj,k,l) 


P 


point 


P 


f(ij+hk,l) 


N 


north 


p+1 


f(ij- hk, I) 


S 


south 


p-1 


f(i+l,j,k,l) 


E 


east 


P + n i 


f{i-hiXl) 


W 


west 


P ~ n j 


f(i,j,k+l,l) 


T 


top 


P + riij 


f(i,j,k-l,t) 


B 


bottom 


p - riij 


f(i,j,k,l+l) 


U 


up 


P + riijk 


f(i,j,k,l-l) 


D 


down 


P ~ n ijk 




Mixed Diagonals 




f(i,j,k + 1,1+1) 


UT 


uptop 


p + (n uk + tlij) 


f(i,j,k- 1,1+1) 


UB 


upbottom 


P + (riijk - fijj) 


f(i,j,k+ 1,1-1) 


DT 


downtop 


P ~ (riijk ~ n ij) 


f(i,j,k- 1,1-1) 


DB 


downbottom 


p-(n ijk + n u ) 
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